function sparse_eigenvalue_test()
% SPARSE_EIGENVALUE_TEST 稀疏矩阵特征值计算专项测试
% 
% 测试大规模稀疏矩阵特征值计算的各种方法和技术

fprintf('=== 稀疏矩阵特征值计算专项测试 ===\n\n');

%% 测试1：稀疏矩阵基础
fprintf('测试1：稀疏矩阵基础\n');
fprintf('目标：理解稀疏矩阵的存储和基本操作\n\n');

% 创建稀疏三对角矩阵
n1 = 100;
e = ones(n1, 1);
A1 = spdiags([e, 2*e, e], [-1, 0, 1], n1, n1);

fprintf('稀疏三对角矩阵: %d×%d\n', n1, n1);
fprintf('非零元素数: %d\n', nnz(A1));
fprintf('稀疏度: %.2f%%\n', 100 * (1 - nnz(A1)/(n1^2)));
fprintf('存储效率: %.2f%% (相对于满矩阵)\n', 100 * nnz(A1)/(n1^2));

% 显示稀疏模式
if n1 <= 20
    fprintf('稀疏模式:\n');
    spy(A1);
    title('稀疏三对角矩阵模式');
end

% 与满矩阵的内存比较
memory_sparse = whos('A1');
A1_full = full(A1);
memory_full = whos('A1_full');

fprintf('内存使用比较:\n');
fprintf('稀疏矩阵: %d bytes\n', memory_sparse.bytes);
fprintf('满矩阵: %d bytes\n', memory_full.bytes);
fprintf('内存节省: %.1f倍\n', memory_full.bytes / memory_sparse.bytes);

%% 测试2：EIGS函数基础
fprintf('\n测试2：EIGS函数基础\n');
fprintf('目标：掌握MATLAB的eigs函数使用\n\n');

% 计算不同数量的特征值
k_values = [1, 3, 6, 10];
fprintf('特征值数量    计算时间    内存使用    精度\n');
fprintf('----------    --------    --------    ----\n');

for k = k_values
    tic;
    [V, D] = eigs(A1, k, 'largestabs');
    time_eigs = toc;
    
    eigenvals = diag(D);
    
    % 验证精度
    residuals = zeros(k, 1);
    for i = 1:k
        residuals(i) = norm(A1 * V(:, i) - eigenvals(i) * V(:, i));
    end
    max_residual = max(residuals);
    
    fprintf('%8d      %8.4f    %8.2f    %.2e\n', ...
            k, time_eigs, memory_sparse.bytes/1e6, max_residual);
end

%% 测试3：不同特征值选择策略
fprintf('\n测试3：不同特征值选择策略\n');
fprintf('目标：比较不同的特征值选择方法\n\n');

k3 = 6;
strategies = {'largestabs', 'smallestabs', 'largestreal', 'smallestreal', 'bothendsreal'};
strategy_names = {'最大模', '最小模', '最大实部', '最小实部', '两端实部'};

fprintf('选择策略        计算时间    特征值范围\n');
fprintf('--------        --------    ----------\n');

for i = 1:length(strategies)
    try
        tic;
        [V3, D3] = eigs(A1, k3, strategies{i});
        time_strategy = toc;
        
        eigenvals3 = diag(D3);
        eigenval_range = [min(real(eigenvals3)), max(real(eigenvals3))];
        
        fprintf('%-12s    %8.4f    [%.4f, %.4f]\n', ...
                strategy_names{i}, time_strategy, eigenval_range);
    catch ME
        fprintf('%-12s    失败: %s\n', strategy_names{i}, ME.message);
    end
end

%% 测试4：位移-逆变换方法
fprintf('\n测试4：位移-逆变换方法\n');
fprintf('目标：使用位移寻找特定区域的特征值\n\n');

% 寻找接近特定值的特征值
target_values = [1.0, 2.0, 3.0];
fprintf('目标值    找到的特征值    误差        计算时间\n');
fprintf('------    ------------    ----        --------\n');

for target = target_values
    tic;
    try
        [V4, D4] = eigs(A1, 3, target);
        time_shift = toc;
        
        eigenvals4 = diag(D4);
        [min_error, min_idx] = min(abs(eigenvals4 - target));
        closest_eigenval = eigenvals4(min_idx);
        
        fprintf('%6.1f    %12.6f    %.2e    %8.4f\n', ...
                target, closest_eigenval, min_error, time_shift);
    catch ME
        fprintf('%6.1f    计算失败: %s\n', target, ME.message);
    end
end

%% 测试5：大规模稀疏矩阵
fprintf('\n测试5：大规模稀疏矩阵\n');
fprintf('目标：测试不同规模稀疏矩阵的计算性能\n\n');

sizes = [500, 1000, 2000, 5000];
fprintf('矩阵规模    非零元素    EIGS时间    EIG时间    加速比\n');
fprintf('--------    --------    --------    -------    ------\n');

for n = sizes
    % 创建稀疏矩阵
    e_large = ones(n, 1);
    A_large = spdiags([e_large, 2*e_large, e_large], [-1, 0, 1], n, n);
    
    % EIGS计算
    k_large = min(10, n-1);
    tic;
    [V_large, D_large] = eigs(A_large, k_large, 'largestabs');
    time_eigs_large = toc;
    
    % 对于小规模矩阵，比较完整特征值分解
    if n <= 1000
        tic;
        eigenvals_full = eig(full(A_large));
        time_eig_full = toc;
        speedup = time_eig_full / time_eigs_large;
    else
        time_eig_full = NaN;
        speedup = NaN;
    end
    
    fprintf('%6d      %8d    %8.4f    %7.4f    %6.1f\n', ...
            n, nnz(A_large), time_eigs_large, time_eig_full, speedup);
end

%% 测试6：不同稀疏结构
fprintf('\n测试6：不同稀疏结构\n');
fprintf('目标：比较不同稀疏结构对计算的影响\n\n');

n6 = 200;
structures = {
    '三对角', @() spdiags([ones(n6,1), 2*ones(n6,1), ones(n6,1)], [-1,0,1], n6, n6);
    '五对角', @() spdiags([ones(n6,1), ones(n6,1), 3*ones(n6,1), ones(n6,1), ones(n6,1)], [-2,-1,0,1,2], n6, n6);
    '随机稀疏', @() sprandsym(n6, 0.05);
    '块对角', @() blkdiag(randn(50,50), randn(50,50), randn(50,50), randn(50,50))
};

fprintf('矩阵结构      非零元素    条件数      EIGS时间\n');
fprintf('--------      --------    ------      --------\n');

for i = 1:size(structures, 1)
    struct_name = structures{i, 1};
    struct_func = structures{i, 2};
    
    try
        A_struct = struct_func();
        if size(A_struct, 1) ~= n6
            A_struct = A_struct(1:min(n6, end), 1:min(n6, end));
        end
        
        % 确保矩阵是对称的
        A_struct = (A_struct + A_struct') / 2;
        
        nnz_count = nnz(A_struct);
        
        % 计算条件数（使用部分特征值估计）
        try
            lambda_max = eigs(A_struct, 1, 'largestreal');
            lambda_min = eigs(A_struct, 1, 'smallestreal');
            cond_est = lambda_max / lambda_min;
        catch
            cond_est = NaN;
        end
        
        % EIGS计算时间
        tic;
        [V_struct, D_struct] = eigs(A_struct, 6, 'largestabs');
        time_struct = toc;
        
        fprintf('%-12s  %8d    %8.2e    %8.4f\n', ...
                struct_name, nnz_count, cond_est, time_struct);
    catch ME
        fprintf('%-12s  计算失败: %s\n', struct_name, ME.message);
    end
end

%% 测试7：收敛性和精度控制
fprintf('\n测试7：收敛性和精度控制\n');
fprintf('目标：分析EIGS的收敛参数对结果的影响\n\n');

A7 = spdiags([ones(200,1), 3*ones(200,1), ones(200,1)], [-1,0,1], 200, 200);

% 不同的容差设置
tolerances = [1e-6, 1e-8, 1e-10, 1e-12];
fprintf('容差        计算时间    最大残差    迭代次数\n');
fprintf('----        --------    --------    --------\n');

for tol = tolerances
    opts.tol = tol;
    opts.disp = 0;  % 关闭显示
    
    tic;
    [V7, D7] = eigs(A7, 5, 'largestabs', opts);
    time_tol = toc;
    
    % 计算残差
    eigenvals7 = diag(D7);
    residuals7 = zeros(5, 1);
    for i = 1:5
        residuals7(i) = norm(A7 * V7(:, i) - eigenvals7(i) * V7(:, i));
    end
    max_residual7 = max(residuals7);
    
    fprintf('%.0e    %8.4f    %.2e    %8s\n', ...
            tol, time_tol, max_residual7, 'N/A');
end

%% 测试8：内存效率分析
fprintf('\n测试8：内存效率分析\n');
fprintf('目标：分析稀疏矩阵计算的内存使用\n\n');

sizes_mem = [100, 500, 1000];
fprintf('矩阵规模    稀疏存储    满矩阵存储    内存节省    计算加速\n');
fprintf('--------    --------    ----------    --------    --------\n');

for n = sizes_mem
    % 创建稀疏矩阵
    A_sparse = spdiags([ones(n,1), 2*ones(n,1), ones(n,1)], [-1,0,1], n, n);
    
    % 内存使用
    info_sparse = whos('A_sparse');
    A_full_temp = full(A_sparse);
    info_full = whos('A_full_temp');
    
    memory_saving = info_full.bytes / info_sparse.bytes;
    
    % 计算时间比较（仅对小规模矩阵）
    if n <= 500
        % 稀疏计算
        tic;
        [V_sp, D_sp] = eigs(A_sparse, 5, 'largestabs');
        time_sparse = toc;
        
        % 满矩阵计算
        tic;
        eigenvals_full_temp = eig(A_full_temp);
        time_full = toc;
        
        compute_speedup = time_full / time_sparse;
    else
        compute_speedup = NaN;
    end
    
    fprintf('%6d      %8.2f    %10.2f    %8.1f    %8.1f\n', ...
            n, info_sparse.bytes/1e6, info_full.bytes/1e6, ...
            memory_saving, compute_speedup);
    
    clear A_full_temp;  % 释放内存
end

%% 测试9：特殊应用：图的拉普拉斯矩阵
fprintf('\n测试9：特殊应用：图的拉普拉斯矩阵\n');
fprintf('目标：计算图拉普拉斯矩阵的特征值\n\n');

% 创建简单图的拉普拉斯矩阵
n_nodes = 50;
% 创建环形图
A_adj = diag(ones(n_nodes-1, 1), 1) + diag(ones(n_nodes-1, 1), -1);
A_adj(1, n_nodes) = 1; A_adj(n_nodes, 1) = 1;  % 闭合环

% 拉普拉斯矩阵 L = D - A
D = diag(sum(A_adj, 2));  % 度矩阵
L = D - A_adj;

fprintf('图节点数: %d\n', n_nodes);
fprintf('图边数: %d\n', nnz(A_adj)/2);
fprintf('拉普拉斯矩阵非零元素: %d\n', nnz(L));

% 计算最小的几个特征值（连通性分析）
[V_lap, D_lap] = eigs(L, 6, 'smallestabs');
eigenvals_lap = diag(D_lap);

fprintf('最小的6个特征值:\n');
for i = 1:6
    fprintf('  λ%d = %.6f\n', i, eigenvals_lap(i));
end

% 第二小特征值（代数连通度）
algebraic_connectivity = eigenvals_lap(2);
fprintf('代数连通度: %.6f\n', algebraic_connectivity);

if algebraic_connectivity > 1e-10
    fprintf('图是连通的\n');
else
    fprintf('图可能不连通\n');
end

%% 测试10：性能基准测试
fprintf('\n测试10：性能基准测试\n');
fprintf('目标：建立稀疏特征值计算的性能基准\n\n');

fprintf('=== 性能基准测试结果 ===\n');
fprintf('测试环境: MATLAB %s\n', version);
fprintf('测试时间: %s\n', datestr(now));

% 标准测试用例
test_cases = {
    '小规模三对角 (100×100)', 100, @(n) spdiags([ones(n,1), 2*ones(n,1), ones(n,1)], [-1,0,1], n, n);
    '中规模三对角 (1000×1000)', 1000, @(n) spdiags([ones(n,1), 2*ones(n,1), ones(n,1)], [-1,0,1], n, n);
    '大规模三对角 (5000×5000)', 5000, @(n) spdiags([ones(n,1), 2*ones(n,1), ones(n,1)], [-1,0,1], n, n);
};

fprintf('\n基准测试结果:\n');
fprintf('测试用例                    矩阵规模    计算时间    特征值精度\n');
fprintf('--------                    --------    --------    ----------\n');

for i = 1:size(test_cases, 1)
    test_name = test_cases{i, 1};
    n_test = test_cases{i, 2};
    matrix_func = test_cases{i, 3};
    
    A_test = matrix_func(n_test);
    k_test = min(10, n_test-1);
    
    tic;
    [V_test, D_test] = eigs(A_test, k_test, 'largestabs');
    time_test = toc;
    
    % 计算精度
    eigenvals_test = diag(D_test);
    residual_test = norm(A_test * V_test(:, 1) - eigenvals_test(1) * V_test(:, 1));
    
    fprintf('%-26s  %8d    %8.4f    %.2e\n', ...
            test_name, n_test, time_test, residual_test);
end

fprintf('\n=== 稀疏矩阵特征值计算测试总结 ===\n');
fprintf('1. 稀疏存储大幅节省内存，特别适合大规模问题\n');
fprintf('2. EIGS函数提供灵活的特征值选择策略\n');
fprintf('3. 位移-逆变换方法可以寻找特定区域的特征值\n');
fprintf('4. 计算时间主要取决于矩阵规模和所需特征值数量\n');
fprintf('5. 不同稀疏结构对计算效率有显著影响\n');
fprintf('6. 容差设置需要平衡精度和计算时间\n');
fprintf('7. 稀疏方法在大规模问题上有巨大优势\n');
fprintf('8. 图论应用中拉普拉斯矩阵特征值有重要意义\n\n');

end